% S-OMP alogrithm for MMV/SMV: Y=AX which comes from the CTF (for real-valued MB signal)
% This function returns a solution Xr together with the support (Supp)
% The termination criteria are:
%    Stop if more than N bands were identified (N is the total number of bands in both parts of the spectrum range)
%    Stop if residual norm is greater than ResThreshold
%    Stop if the ratio between residual norm and the norm of the current solution is greater than ResvsSolThreshold
% Each iteration picks one support item. The symmetric location is added
% automatically. Once the support contains 2N locations, we check  adjacent
% locations. Finally we repeat "MoreIters" (with symmetric completion as
% before).
function [Supp] = SOMP(S,A,N)
    [mA,nA] = size(A);
    [mY,nY] = size(S);

    %Initialize
    residual = S;
%     OneSidedSupp = [];
    Supp = [];
%     SymmetricSupp = [];
    iter = 1;
    NormACols = sqrt(diag(A'*A));

    while (iter <= N)
            % Matching step - solve optimization problem
            Z_1 = A'*residual;
            Z = sqrt(sum(abs(Z_1).^2,2))./NormACols;
            [~, maxPos] = max(Z);
            BestLoc = maxPos(1);

            % Update support - Because of the SBR2 Method the support
            % indecs are not syymetric
            Supp = [Supp BestLoc];
%             SymmetricLoc =  (nA+1-BestLoc);
%             if (BestLoc ~= SymmetricLoc)
%                 SymmetricSupp = [SymmetricSupp SymmetricLoc];
%             end
%             Supp = [OneSidedSupp  SymmetricSupp];

            % Project residual
            As = A(:,Supp);
            solution = As*pinv(As)*S;
            residual = S-solution;

            iter = iter+1;
            
            
    end


end

